************************************************************************************************
* This replication file estimates a free entry counterfactual with a potential price reduction *
************************************************************************************************
*1. Re-calibrate parameters to reflect price reduction.
*2. Estimate and store information on the status quo for comparison.
*3. Estimate the counterfactual distribution of firms under free entry.
*4. Store the outcomes
***********************
* 0 Set path to files *
***********************
clear all
capture log close 
set more off
set trace off
mata: mata set matastrict off
cd $main_directory

**************************************************************************************
* 1. Load the parameters and re-calibrate the intercept to reflect price reductions. *
**************************************************************************************
sca counter_offices=1 //Set this variable to 1 if counterfactual office locations in the addresses of postoffices

/*
// Option 1: Parallel computing
//The variable PBS_ARRAY_INDEX refers to a row index of a matrix, containing the district identifier in the first column, the price reduction for real-estate in the second column and the price reduction for other transactions in the third.
//Thus, the matrix contains all the combinations of district and price reduction that we look at in order to find the optimal price reduction under free entry.
global selection_counterfactual : env PBS_ARRAY_INDEX
mata: mata matuse "Estimates/counterfactual_grid_free", replace
mata:    st_numscalar("district",counterfactual_grid[`=$selection_counterfactual',1])
mata:    st_numscalar("reduction_P1",counterfactual_grid[`=$selection_counterfactual',2])
mata:    st_numscalar("reduction_P2",counterfactual_grid[`=$selection_counterfactual',3])
*/

//Option 2: Manual input of district identifier and price reduction levels
sca district=5 //The 5-th district is Leuven, a relatively small district which computes fairly quickly and can thus be used to test the code.
sca reduction_P1=19 //Price reduction for real estate transactions
sca reduction_P2=19 //Price reduction for other transactions

//Set predetermined values:
global select_arrond `=district'
global reduction_level_Q1 `=reduction_P1'
global reduction_level_Q2 `=reduction_P2'

global starting_value `" 1 "' //If this is set to 1, then we begin with the status quo firm distribution starting values.

use "Data_Belgium\Generated_data\baseline_demand_dataset.dta", clear

if counter_offices==1 {
append using "Data_Belgium\Generated_data\baseline_dataset_counter_post_offices.dta"
}
qui do "Command_files\spatial_demand_Q.do"

//Set predetermined values:


sca scale=10
set seed 1234
mata: mata matuse "Estimates\betas_Q", replace
mata: mata matuse "Estimates\betas_Q1", replace
mata: mata matuse "Estimates\betas_Q2", replace
mata: mata matuse "Estimates\cost_price_parameters", replace
mata:    st_numscalar("p_Q1",p_Q1)
mata:    st_numscalar("p_Q2",p_Q2)
mata:    st_numscalar("alpha",alpha)
mata:    st_numscalar("nu",nu)
mata:    st_numscalar("aL",aL)
mata:    st_numscalar("phi",phi)
mata:    st_numscalar("aM",aM)
mata:    st_numscalar("fixed_cost",fixed_cost)
mata: lambda=1 //This is the weight put on consumer surplus when calculating welfare.


if $select_arrond == 3 {//The third and seventh district are joined together legally.
global select_arrond_secondary `" 7 "' 
}
else{
global select_arrond_secondary `" 0 "' 
}

mata: change_p_Q1=1-`=$reduction_level_Q1'/100
mata: change_p_Q2=1-`=$reduction_level_Q2'/100
//If a price change is planned, make sure that the price change is reflected in the variables 
mata: p_Q1_new=change_p_Q1*p_Q1 
mata: p_Q2_new=change_p_Q2*p_Q2 
//The new estimate vector will be adjusted to reflect lower prices
mata: betas_Q1_new=betas_Q1
mata: betas_Q2_new=betas_Q2 
mata:    st_numscalar("p_Q1_new",p_Q1_new)
mata:    st_numscalar("p_Q2_new",p_Q2_new)
//beta_star reflects how much a Euro of transportation costs is worth in terms of utility
mata: beta_star_Q1=betas_Q1[1]/transportation_cost 
mata: beta_star_Q2=betas_Q2[1]/transportation_cost 
//gamma_star reflects the intercept net of price effects
mata: gamma_star_Q1=betas_Q1[rows(betas_Q1)-1]-beta_star_Q1*p_Q1
mata: gamma_star_Q2=betas_Q2[rows(betas_Q2)-1]-beta_star_Q2*p_Q2
//The new intercept is gamma_star plus the new price effects
mata: betas_Q1_new[rows(betas_Q1_new)-1]=gamma_star_Q1+beta_star_Q1*p_Q1_new
mata: betas_Q2_new[rows(betas_Q2_new)-1]=gamma_star_Q2+beta_star_Q2*p_Q2_new

***********************************************************************
* 2. Estimate and store information on the status quo for comparison. *
***********************************************************************
//Calculate the total demand for real estate on the market
qui su Q1 if common==1
sca sumQ1=r(sum)
sca di sumQ1
//Calculate the total demand for other transactions on the market
qui su Q2 if common==1
sca sumQ2=r(sum)
sca di sumQ2
//Calculate the total population on the market
bysort market_ID:  gen dup_market = cond(_N==1,0,_n)
qui su population if dup_market<2
sca sumPopulation=r(sum)
sca di sumPopulation
//Calculate total potential demand per capita as 10 times larger than the current value.
scalar gamma_Q1=scale*sumQ1/sumPopulation
sca di gamma_Q1
scalar gamma_Q2=scale*sumQ2/sumPopulation
sca di gamma_Q2
//Calculate the market size (i.e. number of possible transactions)
gen market_size_Q1=gamma_Q1*population
gen market_size_Q2=gamma_Q2*population

//Define the product for which demand will be estimated.
gen firm_demand_Q1=Q1 
gen firm_demand_Q2=Q2

//Generate a unique numeric identifier for the market.
sort market_ID
egen market_ID_num = group(market_ID)
//Keep only observations in the selected judicial district
keep if arrondjud==$select_arrond | arrondjud==$select_arrond_secondary

// Controls
global controls_utility distance lnNotPerOff start_not LDE Perc_Young Perc_Elderly log_income popdens
global controls_total_Q1 newid firm_demand_Q1 market_ID_num market_size_Q1 $controls_utility
global controls_total_Q2 newid firm_demand_Q2 market_ID_num market_size_Q2 $controls_utility

//Define newid, which is just firm_IDs which go from 1 to 1152.
egen newid = group(firm_ID)
// Generate a unique observation identifier, which is used to merge the results from the counterfactuals with the initial dataset
tostring market_ID_num, gen(market_ID_string)
tostring newid, gen(newid_string)
gen join="_"
egen unique_identifier=concat(newid_string join market_ID_string)
drop join market_ID_string newid_string

// Count the number of firms in the sample
bysort firm_ID:  gen dup_firm = cond(_N==1,0,_n)
count if dup_firm<2
sca n_firm=r(N)
// Count the number of markets in the subsample
bysort market_ID: gen nvals = _n == 1 
count if nvals
sca markets_number=r(N)
di markets_number
drop nvals
* Count the number of markets with a local notary
preserve
keep if common==1 //This denotes the local market of a notary
bysort market_ID: gen nvals = _n == 1 
count if nvals
sca covered_markets_current=r(N)
sca uncovered_markets_current=markets_number-covered_markets_current
su NotPerOff
sca total_notaries_current=r(sum)
restore

//STATUS QUO ESTIMATION AND DESCRIPTIVES
* Sort the data by market ID, since this will be the first relevant grouping (to calculate notary utility within a municipality)
sort market_ID_num newid
mata:
// Select the relevant variables and read them into mata
X_Q1    = st_data(., "$controls_total_Q1")
X_Q2    = st_data(., "$controls_total_Q2")
// nx counts the number of observations
nx    = rows(X_Q1)
// Add an intercept, with length equal to the number of observations
X_Q1    = X_Q1,J(nx,1,1)
X_Q2    = X_Q2,J(nx,1,1)

// Calculate per firm profit as a matrix 
n_firm=st_numscalar("n_firm")
total_notaries_current=st_numscalar("total_notaries_current")
Results_ln_current_Q1=spdemand_log(X_Q1,betas_Q1)
X_Q1=sort(X_Q1,(3,1)) //sort by market and then firm_id
Results_ln_current_Q2=spdemand_log(X_Q2,betas_Q2)
X_Q2=sort(X_Q2,(3,1)) //sort by market and then firm_id
Results_current_Q1=exp(Results_ln_current_Q1[,2])
Results_current_Q2=exp(Results_ln_current_Q2[,2])
Results_current=Results_current_Q1+Results_current_Q2
output_current=round(sum(Results_current))
consumer_surplus_market_current=consumer_surplus(X_Q1,betas_Q1,transportation_cost)+consumer_surplus(X_Q2,betas_Q2,transportation_cost)
per_firm_variable_profit_current=per_firm_profit_mult_v2(X_Q1,X_Q2,betas_Q1,betas_Q2,p_Q1,p_Q2,nu,aL,phi,aM,alpha)
consumer_surplus_current=round(sum(consumer_surplus_market_current[,2]))
producer_surplus_current=round(sum(per_firm_variable_profit_current[,2])-fixed_cost*total_notaries_current)
welfare_current=round(consumer_surplus_current+producer_surplus_current)
   st_numscalar("output_current", output_current)
   st_numscalar("consumer_surplus_current", consumer_surplus_current)
   st_numscalar("producer_surplus_current", producer_surplus_current)
   st_numscalar("welfare_current", welfare_current)
//Generate a vector with the number of notaries in each office for every office-market combination
n_notaries=J(nx,1,0)
n_notaries=round(exp(X_Q1[.,6])) //The sixth column of X stores the log of the number of notaries
//Generate a vector with the number of notaries for every office
notary_vector=J(n_firm,1,0) //Make an empty notary vector
//Calculate the actual notaries per office under the status quo
for(f=1; f<=n_firm; f++) {
notary_vector[f,1]=sum(n_notaries:*(X_Q1[.,1]:==f))/(sum((X_Q1[.,1]:==f))) //This is essentially calculating the average number of notaries across all observations for the same notary office, which is just the number of notaries for that office.
}

//Update the current variable profit to the reduced price
per_firm_variable_profit_current=per_firm_profit_mult_v2(X_Q1,X_Q2,betas_Q1_new,betas_Q2_new,p_Q1_new,p_Q2_new,nu,aL,phi,aM,alpha)
firm_profit_old=per_firm_variable_profit_current[.,2]-fixed_cost*notary_vector //This is the vector of current profits, which will be updated with each round of added/removed notaries.
firm_profit_old[.,1]=mm_cond((firm_profit_old[.,1]:==.),0,firm_profit_old[.,1]) //If firm profits are missing, this means that there are no notaries in the office, so we replace firm profits with 0 in that case.
firm_profit_new=J(n_firm,n_firm,0) //This is a vector holding the counterfactual information on firm-level profits if entry occurs in a given location.
overall_change=100 //This variable is going to measure whether an extra loop of calculations of profits is necessary to check for the need of additional entry or exit. We set it to a positive number (e.g. 100) at the beginning, in order to make sure that we perform at least one round of checking for potential entry and exit under the new counterfactual.
round=0
	end
di consumer_surplus_current
di producer_surplus_current
di welfare_current

*****************************************
* 3. Run the counterfactual estimation. *
*****************************************
//0. Store the firm distribution at the start of the loop. We will compare this to the notary distribution at the end of the loop to know if we should keep running the counterfactual.
//1. Checking need for removal:
//1.1 Calculate per-firm profits conditional on the distribution of firms at the start of the loop and store this as firm_profit_old and per_notary_profit_old, since these are the profits relative to which we will be calculating returns from removing a notary.
//1.2 If per-notary profits are negative, remove a notary from the least-profitable location (based on per-notary profits stored in per_notary_profit_old).
//1.3 Calculate per-firm and per-notary profits after the removal of the least profitable notary. If some firms are still unprofitable, keep repeating step 1.2 and 1.3 until all unprofitable notaries have been removed.
//2. Checking the need for addition:
//2.1 Add a notary to each possible location and calculate the expected profits after the addition. Store the profits of all notaries into a matrix called firm_profit_new. The diagonal of this matrix are the profits in the office where the notary was introduced.
//2.2 Calculte per-notary profits at each office after a notary has been added to it. Find the most profitable location, check that addition increases the firm's profits and add a notary, if that is the case.
mata: 
while (overall_change!=0) { //This checks if after entry there was also exit, meaning that we need to check for positive returns to entry again.
//0. Store the firm distribution at the start of the loop. We will compare this to the notary distribution at the end of the loop to know if we should keep running the counterfactual.
notary_vector_loop_start=notary_vector
//****************************************
//* STEP 1: REMOVE UNPROFITABLE NOTARIES *
//****************************************
//1.1 Calculate per-firm profits conditional on the distribution of firms at the start of the loop and store this as firm_profit_old and per_notary_profit_old, since these are the profits relative to which we will be calculating returns from removing.
//Select observations with a non-zero number of notaries
X_selection1_Q1=select(X_Q1, X_Q1[.,6]:>=0) //Select those values of X where the log of notaries is >=0, i.e. only use data on active notary offices.
X_selection_Q1=select(X_selection1_Q1, X_selection1_Q1[.,6]:!=.) //Select those values of X where the log of notaries is not missing, i.e. only use data on active notary offices.
X_selection1_Q2=select(X_Q2, X_Q2[.,6]:>=0) //Select those values of X where the log of notaries is >=0, i.e. only use data on active notary offices.
X_selection_Q2=select(X_selection1_Q2, X_selection1_Q2[.,6]:!=.) //Select those values of X where the log of notaries is not missing, i.e. only use data on active notary offices.

//Calculate variable profits per firm
per_firm_variable_profit=per_firm_profit_mult_v2(X_selection_Q1,X_selection_Q2,betas_Q1_new,betas_Q2_new,p_Q1_new,p_Q2_new,nu,aL,phi,aM,alpha) //Calculate per firm variable profits

//Make an matrix in which  to store the profits of the firms. The first column contains the id of the firm, while the second column contains the profits (initially set to zero and then filled in step by step for those notaries who are active).
//Make the matrices which are to be filled in. The first column contains the id of the firm, while the second column contains the profits (initially set to zero).
per_firm_profit=range(1,rows(notary_vector),1),J(rows(notary_vector),1,0)
per_notary_profit=range(1,rows(notary_vector),1),J(rows(notary_vector),1,0)
//Fill in the profits of active firms 
for (f=1; f<=rows(per_firm_variable_profit); f++){ //f goes over all active notaries
id_of_firm=per_firm_variable_profit[f,1]
//Calculate firm profits. Note that the notary vector and per-firm variable profits are assumed to be ordered by the id of the firm.
per_firm_profit[id_of_firm,2]=per_firm_variable_profit[f,2]-fixed_cost*notary_vector[id_of_firm,1] //Subtract the fixed costs to obtain total firm profits
per_notary_profit[id_of_firm,2]=(per_firm_variable_profit[f,2]-fixed_cost*notary_vector[id_of_firm,1])/notary_vector[id_of_firm,1] //Calculate per notary profits
} //End of loop filling in available profits
//Store firm profits as firm_profit_old.
firm_profit_old=per_firm_profit[.,2]
per_notary_profit_old=per_notary_profit[.,2]

//1.2 If per-notary profits are negative, remove a notary from the least-profitable location (based on per-notary profits).
//Beginning of loop which checks if removing a notary from a location raises profits.
while (sum(per_notary_profit_old:<0)>0) { //The loop goes on, as long as some offices make negative profits.
minindex(per_notary_profit_old, 1, min_ID=., junk=.) //This returns the ID of the least profitable location.
min_ID=min_ID[1,1] //In case of multiple offices whose removal increases profits by the same amount, we remove the one which is ordered first.
notary_vector[min_ID,1]=notary_vector[min_ID,1]-1 //Remove a notary from the unprofitable location.
X_Q1[.,6]=mm_cond((X_Q1[.,1]:==min_ID),ln(notary_vector[min_ID,1]),X_Q1[.,6]) //Change the control variable matrix to reflect this change.
X_Q2[.,6]=mm_cond((X_Q2[.,1]:==min_ID),ln(notary_vector[min_ID,1]),X_Q2[.,6]) //Change the control variable matrix to reflect this change.

//Step 1.3: Calculate profits after the removal of the notary, i.e. calculation of profits for all offices when the chosen office has n-1 notaries
//Select observations with a non-zero number of notaries
X_selection1_Q1=select(X_Q1, X_Q1[.,6]:>=0) //Select those values of X where the log of notaries is >=0, i.e. only use data on active notary offices.
X_selection_Q1=select(X_selection1_Q1, X_selection1_Q1[.,6]:!=.) //Select those values of X where the log of notaries is not missing, i.e. only use data on active notary offices.
X_selection1_Q2=select(X_Q2, X_Q2[.,6]:>=0) //Select those values of X where the log of notaries is >=0, i.e. only use data on active notary offices.
X_selection_Q2=select(X_selection1_Q2, X_selection1_Q2[.,6]:!=.) //Select those values of X where the log of notaries is not missing, i.e. only use data on active notary offices.

//Calculate variable profits per firm
per_firm_variable_profit=per_firm_profit_mult_v2(X_selection_Q1,X_selection_Q2,betas_Q1_new,betas_Q2_new,p_Q1_new,p_Q2_new,nu,aL,phi,aM,alpha) //Calculate per firm variable profits
//Make the matrices which are to be filled in. The first column contains the id of the firm, while the second column contains the profits (initially set to zero).
per_firm_profit=range(1,rows(notary_vector),1),J(rows(notary_vector),1,0)
per_notary_profit=range(1,rows(notary_vector),1),J(rows(notary_vector),1,0)
//Fill in the profits of active firms 
for (f=1; f<=rows(per_firm_variable_profit); f++){ //f goes over all active notaries
id_of_firm=per_firm_variable_profit[f,1]
//Calculate firm profits. Note that the notary vector and per-firm variable profits are assumed to be ordered by the id of the firm.
per_firm_profit[id_of_firm,2]=per_firm_variable_profit[f,2]-fixed_cost*notary_vector[id_of_firm,1] //Subtract the fixed costs to obtain total firm profits
per_notary_profit[id_of_firm,2]=(per_firm_variable_profit[f,2]-fixed_cost*notary_vector[id_of_firm,1])/notary_vector[id_of_firm,1] //Calculate per-notary profits
} //End of loop filling in available profits
//Store firm profits as firm_profit_old. In case some profits are still negative. If that is the case, the loop will go over all offices again (Step 1.2 and 1.3) and continue enacting closures.
firm_profit_old=per_firm_profit[.,2]
per_notary_profit_old=per_notary_profit[.,2]
} //End of if loop checking if notary removal is necessary.

//***********************************
//* STEP 2: ADD PROFITABLE NOTARIES *
//***********************************
//2.1 Add a notary to each possible location and calculate the expected profits after the addition.
extra_profit=100 //This variable denotes the profit of the last entrant. The value of 100 will be replaced with actual profits during the first round of the loop. It simply ensures that we check at least once for a profit-increasing entry.
//Beginning of loop which checks if adding a notary to a location is profitable
while (extra_profit>0) { //The loop goes on, as long as entry raises the profits of the office where the entry occurs.
for(i=1; i<=n_firm; i++) { //i goes over the notaries, adds a notary and calculates profits after addition.
//Add a notary to office i
notary_vector[i,1]=	notary_vector[i,1]+1 //Add the notary to the notary vector where the i-th element shows how many notaries are active in location i.
X_Q1[.,6]=mm_cond((X_Q1[.,1]:==i),ln(notary_vector[i,1]),X_Q1[.,6]) //Change the control variable matrix to reflect this change.
X_Q2[.,6]=mm_cond((X_Q2[.,1]:==i),ln(notary_vector[i,1]),X_Q2[.,6]) //Change the control variable matrix to reflect this change.

//Calculation of profits for all offices when the chosen office has n+1 notaries
//Select observations with a non-zero number of notaries
X_selection1_Q1=select(X_Q1, X_Q1[.,6]:>=0) //Select those values of X where the log of notaries is >=0, i.e. only use data on active notary offices.
X_selection_Q1=select(X_selection1_Q1, X_selection1_Q1[.,6]:!=.) //Select those values of X where the log of notaries is not missing, i.e. only use data on active notary offices.
X_selection1_Q2=select(X_Q2, X_Q2[.,6]:>=0) //Select those values of X where the log of notaries is >=0, i.e. only use data on active notary offices.
X_selection_Q2=select(X_selection1_Q2, X_selection1_Q2[.,6]:!=.) //Select those values of X where the log of notaries is not missing, i.e. only use data on active notary offices.
//Calculate per-firm profits after the addition.
per_firm_variable_profit=per_firm_profit_mult_v2(X_selection_Q1,X_selection_Q2,betas_Q1_new,betas_Q2_new,p_Q1_new,p_Q2_new,nu,aL,phi,aM,alpha) //Calculate per firm variable profits

//Make an matrix in which  to store the profits of the firms. The first column contains the id of the firm, while the second column contains the profits (initially set to zero).
per_firm_profit=range(1,rows(notary_vector),1),J(rows(notary_vector),1,0)
//Fill in the profits of active firms 
for (f=1; f<=rows(per_firm_variable_profit); f++){
id_of_firm=per_firm_variable_profit[f,1]
per_firm_profit[id_of_firm,2]=per_firm_variable_profit[f,2]-fixed_cost*notary_vector[id_of_firm,1] //Subtract the fixed costs to obtain total firm profits
} //End of loop filling in available profits
//Store firm profits as a row in firm_profit_new, this will contain all the relevant information necessary to decide where to add a notary. An element i,j of the matrix contains the profits of the j-th notary office if a notary is added to firm i.
firm_profit_new[i,.]=per_firm_profit[.,2]'
//Remove the added notary from all i observations in order to get back to the original notary vector
notary_vector[i,1]=	notary_vector[i,1]-1 //Remove the notary to the notary vector where the i-th element shows how many notaries are active in location i.
X_Q1[.,6]=mm_cond((X_Q1[.,1]:==i),ln(notary_vector[i,1]),X_Q1[.,6]) //Change the control variable matrix to reflect this change.
X_Q2[.,6]=mm_cond((X_Q2[.,1]:==i),ln(notary_vector[i,1]),X_Q2[.,6]) //Change the control variable matrix to reflect this change.
} //End of loop which goes over all offices, adds 1 notary and calculates expected profits.

//2.2 Calculte per-notary profits at each office after a notary has been added to it. Find the most profitable location, check that addition increases the firm's profits and add a notary, if that is the case.
//The diagonal of firm_profit_new contains the profits of notary offices if a notary is added to them.
profit_entry=diagonal(firm_profit_new):/(notary_vector+J(rows(notary_vector),1,1)) //Calculate the per-notary profits after the addition.
//Find the location at which adding a notary leads to the largest increase in profits.
maxindex(profit_entry, 1, max_ID=., junk=.) //This returns the ID of the location where profits are highest per notary after the addition.
max_ID=max_ID[1,1]
//Calculate the extra profit generated from the addition to the most profitable notary office.
extra_profit=profit_entry[max_ID,1]
//Add one notary from the max_ID location, as long as it increases profits.
if (extra_profit>0) { //Start of condition checking that profits from entry are positive and adding a notary to the most profitable location.
notary_vector[max_ID,1]=notary_vector[max_ID,1]+1 //Add the notary to the notary vector where the i-th element shows how many notaries are active in location i.
X_Q1[.,6]=mm_cond((X_Q1[.,1]:==max_ID),ln(notary_vector[max_ID,1]),X_Q1[.,6]) //Change the control variable matrix to reflect this change.
X_Q2[.,6]=mm_cond((X_Q2[.,1]:==max_ID),ln(notary_vector[max_ID,1]),X_Q2[.,6]) //Change the control variable matrix to reflect this change.
firm_profit_old=firm_profit_new[max_ID,.]' //Store the expected profits if entry happens at max_ID as the new current profits.
} //End of condition checking that profits from entry are positive and adding a notary to the most profitable location.
} //End of loop checking if there is any need for entry

//*****************************
//* STEP 4: CHECK CONVERGENCE *
//*****************************
//If there was change in the number of notaries, we need to check again if exit or entry is profitable (i.e. checking Nash-equilibrium conditions).
notary_vector_change=notary_vector-notary_vector_loop_start
overall_change=sum(notary_vector_change:!=0)
}

newid=X_Q1[.,1]
market_ID_num=X_Q1[.,3]
n_notaries=exp(X_Q1[.,6])
join=J(rows(X_Q1),1,"_")
unique_identifier=strofreal(newid)+join+strofreal(market_ID_num)
	end
	

************************************************
* 4. Store the outcomes of the counterfactual. *
************************************************

getmata n_notaries, id(unique_identifier)
mata: notary_vector
mata: zero_notaries=sum((notary_vector[.,1]:==0))
mata: zero_notaries
mata: total_notaries=sum(notary_vector)
mata: total_notaries
mata: offices_number=sum((notary_vector[.,1]:>0))
su NotPerOff
replace NotPerOff=n_notaries
su NotPerOff
su n_notaries
//Keep only active notary offices
drop if NotPerOff==0
drop if NotPerOff==.
replace lnNotPerOff=ln(NotPerOff)

//Calculate consumer surplus and producer surplus
sort market_ID newid
mata:
// Select the relevant variables and read them into mata
X_Q1    = st_data(., "$controls_total_Q1")
X_Q2    = st_data(., "$controls_total_Q2")
// nx counts the number of observations
nx    = rows(X_Q1)
// Add an intercept, with length equal to the number of observations and 
X_Q1    = X_Q1,J(nx,1,1)
X_Q2    = X_Q2,J(nx,1,1)

// Calculate per firm profit as a matrix 
Results_ln_Q2=spdemand_log(X_Q2,betas_Q2_new)
Results_ln_Q1=spdemand_log(X_Q1,betas_Q1_new)
Results_Q1=exp(Results_ln_Q1[,2])
Results_Q2=exp(Results_ln_Q2[,2])
Results=Results_Q1+Results_Q2
mean(Results_Q1)
mean(Results_Q2)
output=round(sum(Results))
consumer_surplus_per_market=consumer_surplus(X_Q1,betas_Q1_new,transportation_cost)+consumer_surplus(X_Q2,betas_Q2_new,transportation_cost)
per_firm_variable_profit=per_firm_profit_mult_v2(X_Q1,X_Q2,betas_Q1_new,betas_Q2_new,p_Q1_new,p_Q2_new,nu,aL,phi,aM,alpha)
consumer_surplus=round(sum(consumer_surplus_per_market[,2]))
producer_surplus=round(sum(per_firm_variable_profit[,2])-fixed_cost*total_notaries)
welfare=round(consumer_surplus+producer_surplus)
st_numscalar("output", output)
st_numscalar("consumer_surplus", consumer_surplus)
st_numscalar("producer_surplus", producer_surplus)
st_numscalar("welfare", welfare)
unprofitable_firms=sum((firm_profit_old:<0))
st_numscalar("unprofitable_firms", unprofitable_firms)
end

keep if common==1 //This denotes the local market of a notary
su NotPerOff
sca total_notaries=r(sum)
sca offices_number=r(N)
bysort market_ID: gen nvals = _n == 1 
count if nvals
sca covered_markets=r(N)
sca uncovered_markets=markets_number-covered_markets
sca change_notaries=total_notaries-total_notaries_current
sca change_offices=offices_number-n_firm
sca change_uncovered_markets=uncovered_markets-uncovered_markets_current
sca change_output=output-output_current
sca change_output_perc=change_output/output_current
sca change_consumer_surplus=consumer_surplus-consumer_surplus_current
sca change_producer_surplus=producer_surplus-producer_surplus_current
sca change_welfare=welfare-welfare_current
//Current outcome
di total_notaries_current " & " n_firm " & " uncovered_markets_current " & " output_current " & " consumer_surplus_current " & " producer_surplus_current " & " welfare_current
//Counterfactual outcome
di total_notaries " & " offices_number " & " uncovered_markets " & " output " & " consumer_surplus " & " producer_surplus " & " welfare
//Changes in outcomes
di change_notaries " & " change_offices " & " change_uncovered_markets " & " change_output " & " change_consumer_surplus " & " change_producer_surplus " & " change_welfare
//Table
di total_notaries " & " offices_number " & " uncovered_markets " &  " output " & " change_consumer_surplus " & " change_producer_surplus " & " change_welfare " & " unprofitable_firms
matrix free_entry=($select_arrond , total_notaries,offices_number,uncovered_markets,output,change_consumer_surplus,change_producer_surplus,change_welfare,unprofitable_firms)

if counter_offices==1 {
mata: free_`=$reduction_level_Q1'_`=$reduction_level_Q2'_`=$select_arrond'_0=st_matrix("free_entry")
//Store the results as a mata table	
mata: mata matsave "Estimates/Entry_counterfactual/free_`=$reduction_level_Q1'_`=$reduction_level_Q2'_`=$select_arrond'_0" free_`=$reduction_level_Q1'_`=$reduction_level_Q2'_`=$select_arrond'_0, replace
//Store the results as a dta file
rename NotPerOff NotPerOff_free_`=$reduction_level_Q1'_`=$reduction_level_Q2'_0
keep bvdid_num NotPerOff_free_`=$reduction_level_Q1'_`=$reduction_level_Q2'_0
save "Estimates/Entry_counterfactual/free_`=$reduction_level_Q1'_`=$reduction_level_Q2'_`=$select_arrond'_0.dta", replace
}
else{
mata: free_`=$reduction_level_Q1'_`=$reduction_level_Q2'_`=$select_arrond'=st_matrix("free_entry")
//Store the results as a mata table	
mata: mata matsave "Estimates/Entry_counterfactual/free_`=$reduction_level_Q1'_`=$reduction_level_Q2'_`=$select_arrond'" free_`=$reduction_level_Q1'_`=$reduction_level_Q2'_`=$select_arrond', replace
//Store the results as a dta file
rename NotPerOff NotPerOff_free_`=$reduction_level_Q1'_`=$reduction_level_Q2'
keep bvdid_num NotPerOff_free_`=$reduction_level_Q1'_`=$reduction_level_Q2'
save "Estimates/Entry_counterfactual/free_`=$reduction_level_Q1'_`=$reduction_level_Q2'_`=$select_arrond'.dta", replace
}
